!
! Kombine     Combine a set of frames.
!
! Copyright © 1998-2006, 2011 F.Hroch (hroch@physics.muni.cz)
!
! This file is part of Munipack.
!
! Munipack is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
! 
! Munipack is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with Munipack.  If not, see <http://www.gnu.org/licenses/>.
!

program skombine

!
! Experimental version, subpixel transformation
!
! Sum (combine) the frames.
!
!

  use interpol
!  use optfile
  use stat
  use robustmean
!  use daoio

  implicit none

! Compiled-in parameters (perhaps, no need to adjust)

! Version:
  Character(len=*), parameter :: version = 'SKOMBINE Version 0.5.x, (C)1997-11 F.Hroch, Masaryk University,Brno,CZ '

! Minimum of values
  Real :: minvalue = 0.0

! filename and line lenght parameters
  integer, parameter :: lenfilename = 1024, lenline = 80  

! max stars to intensity computing
  integer, parameter :: maxstar = 256

! FITS
  Logical :: simple,extend,anyf,exist
  Integer :: n, nalldate, maxalldate
  Integer :: i, j, ii, istat, bpix, maxvalue, bitpix = -32, naxis, naxes(2) 
  Integer :: BlockSize,pcount,gcount
  character(len=8) :: FILTER='FILTER',DATE_OBS='DATE-OBS',EXPTIME='EXPTIME', &
       OBJECT='OBJECT'

  Real, dimension(:,:), allocatable :: KOMBINE, IMAGE
  real, dimension(:), allocatable :: obr
  Real :: bscale, bzero
  Real :: oneexp, etime
  Real :: xx,yy,c,r(2),m(2,2),t1(2),t2(2),x,y,f,fs,sky,amag
  real, dimension(maxstar) :: fx, fx0, asky, asky0
  integer :: ix, iy, nl,eq,cc
  character(len=lenline) :: onedate, obj, filtr
  character(len=lenline), dimension(:), allocatable :: alldate, alld
  character(len=lenfilename) :: sb1,sb2,sb,line,record!,sbtra,sbmat
  character(len=lenline) :: buf,key,value
  logical :: initialize = .true.

! selection of method of interpolation
! 0=no interpolation, 1=bilinear, 2=bicubic, 3=bspline
!   fast                slow        slower     hung up
!   tooths              good        perfect    little bit over perfect
  integer :: sint = 0

! set rescale factor (rescale*dimensions of first image, if nx,ny < 0)
  real :: rescale = 1

! set dimensions of output image 
  Integer ::  nx = -1, ny = -1
  
! Output image name:
  Character(len=lenfilename) :: namekombi='skombine.fits'

! Flux recompute to one level
  logical :: flux = .false.

! Normalize output,
! 0.. no normalization, 1 .. mean (sum divided by number), 2.. normlize to sky
  integer :: norm = 0

! Time start of first exposure
  character(len=lenline) :: first = ''

!------

!  write(*,*) trim(version)

! init variables
  istat = 0
  obj = ''
  etime = 0.0
  nalldate = 0
  maxalldate = 2
  f = 1.0
  fs = 0.0
  sky = 0.0
  allocate(alldate(maxalldate))

! check options file
!!$  if( opt_open(21,'kombine.opt',i) )then
!!$     do
!!$        call opt_read(21,key,value,i)
!!$        if( i == opt_end .or. i == opt_err ) exit
!!$        !write(*,*) trim(key),':',trim(value)
!!$
!!$        if( key == 'Method of interpolation' ) then
!!$           read(value,*,iostat=i) sint
!!$        endif
!!$
!!$        if( key == 'Rescale factor' ) then
!!$           read(value,*,iostat=i) rescale
!!$        endif
!!$        
!!$        if( key == 'Output normalized by' ) then
!!$           read(value,*,iostat=i) norm
!!$        endif
!!$
!!$        if( key == 'Multiply flux' ) then
!!$           read(value,*,iostat=i) flux
!!$        endif
!!$        
!!$        if( key == 'Output width' ) then
!!$           read(value,*,iostat=i) nx
!!$        endif
!!$
!!$        if( key == 'Output height' ) then
!!$           read(value,*,iostat=i) ny
!!$        endif
!!$
!!$        if( key == 'Output BITPIX' ) then
!!$           read(value,*,iostat=i) bitpix
!!$        endif
!!$
!!$        if( key == 'Output filename' ) then
!!$           namekombi = value
!!$        endif
!!$
!!$        if( key == 'Key DATE-OBS' ) then
!!$           read(value,*,iostat=i) DATE_OBS
!!$        endif
!!$
!!$        if( key == 'Key EXPTIME' ) then
!!$           read(value,*,iostat=i) EXPTIME
!!$        endif
!!$
!!$        if( key == 'Key FILTER' ) then
!!$           read(value,*,iostat=i) FILTER
!!$        endif
!!$
!!$        if( key == 'Key OBJECT' ) then
!!$           read(value,*,iostat=i) OBJECT
!!$        endif
!!$
!!$     enddo
!!$     call opt_close(21)
!!$  end if

  ! input of top line with init values
!  read(*,*,end=666,iostat=j) sb,x,bpix,n,i,ii
!  if( j /= 0 ) stop 'Bad format of input data.'
!  if( sb /= '' ) namekombi = sb
!  if( x > epsilon(0.0) ) rescale = x
!  if( bpix /= 0 ) bitpix = bpix
!  if( n >= 0 ) norm = n
!  if( i >= 0 ) sint = i
!  if( ii > 0 ) flux = .true.
  
  n = 0
  istat = 0
! input filename
  Do
     read(*,'(a)',end=20,iostat=i) record
     if( i /= 0 ) stop 'Bad format of input filename.'
     write(*,*) trim(record),':'

     cc = index(record,'#')
     if( cc == 0 ) then
        line = record
     else
        line = record(1:cc-1)
     end if

     eq = index(line,'=')

     if( index(line,'OUTPUT') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) namekombi

     else if( index(line,'RESCALE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) rescale

     else if( index(line,'INTERPOL') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) sint

     else if( index(line,'FLUX') > 0 .and. eq > 0) then

        flux = .true.

     else if( index(line,'FILE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) sb

         
! input initialisation of the current file
     Call ftopen(25,sb,0,BlockSize,istat)
     Call ftghpr(25,2,simple,bpix,naxis,naxes,pcount,gcount,extend,istat)
     if( istat /= 0 ) goto 666
     CALL ftgkye(25,'BSCALE',bscale,buf,istat)
     CALL ftgkye(25,'BZERO',bzero,buf,istat)
     !        BSCALE, BZERO not presented
     if( istat /= 0 ) then
        istat = 0
        bscale = 1.0
        bzero = 0.0
     endif
     call ftgkye(25,EXPTIME,oneexp, buf, istat)
     if( istat /= 0 ) then
        oneexp = -1.0
        istat = 0
     endif
     call ftgkys(25,DATE_OBS,onedate, buf, istat)
     if( istat /= 0 ) then
        onedate = ''
        istat = 0
     endif
     call ftgkys(25,OBJECT,obj, buf, istat)
     if( istat /= 0 ) then
        obj = ''
        istat = 0
     endif
     call ftgkys(25,FILTER,filtr, buf, istat)
     if( istat /= 0 ) then
        filtr = ''
        istat = 0
     endif

! image - data input
     allocate(IMAGE(naxes(1),naxes(2)))
     !IMAGE = 0.0
     call FTG2DE(25,1,minvalue,naxes(1),naxes(1),naxes(2),IMAGE,anyf,istat)
     if( istat /= 0 ) goto 666
     call ftclos(25,istat)
     if( istat /= 0 ) goto 666

! image - transformation input
!!$     sbtra = sb
!!$     call another(sbtra,'.fits','.TRA')
!!$     open(27,file=sbtra,status='old',iostat=i)
!!$     if( i == 0 ) then
!!$        read(27,*,iostat=i) sb1,sb2,c,r,m,t1,t2
!!$        close(27)
!!$     endif
!!$     if( i /= 0 ) then
!!$        write(*,*) 'No such TRA file or invalid format. Skipped.'
!!$        goto 665
!!$     endif
     c = 1.0
     r = 0.0
     m(1,:) = (/1.0,0.0/)
     m(2,:) = (/0.0,1.0/)
     t1 = 0.0
     t2 = 0.0

! image - determine flux parameters
!!$     if( flux ) then
!!$        
!!$        sbmat = sb
!!$        call another(sbmat,'.fits','.MAT')
!!$        j = 0
!!$        call dao_rdhead(28,sbmat,nl,j)
!!$        i = 0
!!$        do
!!$           call dao_rdstar(28, nl, ii, x, y, amag, xx, j)
!!$           if( j == dao_end ) exit
!!$           !write(*,*) i, x, y, amag(i), asky(i)
!!$           if( j == dao_ok ) then
!!$              i = i + 1
!!$              fx(i) = 10**(-0.4*(amag - 25.0))
!!$              asky(i) = xx
!!$           endif
!!$           if( i == maxstar ) exit
!!$        enddo
!!$        call dao_close(28,j)
!!$        if( i > 0 .and. .not. initialize ) then
!!$           where( fx > 0 )
!!$              fx = fx0/fx
!!$           end where
!!$           
!!$           call rmean(fx(1:i),f,fs)
!!$           asky = f*asky - asky0
!!$           call rmean(asky(1:i),sky,x)
!!$           !write(*,*) f,sky
!!$        else
!!$           f = 1.0
!!$           fs = 0.0
!!$           sky = 0.0
!!$        end if
!!$     end if

! first image - internal buffers initialisation
     if( initialize ) then
        initialize = .false.

        if( nx < 0 .and. ny < 0 )then
           nx = nint(rescale*naxes(1))
           ny = nint(rescale*naxes(2))
        endif

        allocate(KOMBINE(nx,ny))
        KOMBINE = 0.0

        if( flux ) then
           fx0 = fx
           asky0 = asky
           f = 1.0
           fs = 0.0
           sky = 0.0
        end if

     end if

     if( sint == 3 ) then
        ! init bspline
        call bspline_init(IMAGE)
     endif

     ! compute sumation
     do i = 1, nx
        do j = 1, ny
           
           xx = i/rescale - t1(1)
           yy = j/rescale - t1(2)
           x = c*(r(1) + xx*m(1,1) + yy*m(1,2)) + t2(1)
           y = c*(r(2) + xx*m(2,1) + yy*m(2,2)) + t2(2)

           !write(*,*) i,j,x,y
              
           select case (sint)

           case(0)
              ! no interpolation
              ix = nint(x)
              iy = nint(y)
              if(  1 <= ix .and. ix <= size(IMAGE,1) .and. &
                   1 <= iy .and. iy <= size(IMAGE,2) )then
                 KOMBINE(i,j) = KOMBINE(i,j) + f*IMAGE(ix,iy)-sky
              end if
                         
           case(1)
              ! bilinear interpolation
              KOMBINE(i,j) = KOMBINE(i,j) + f*bilinear(x,y,IMAGE)-sky
              
           case(2)
              ! bicubic interpolation
              KOMBINE(i,j) = KOMBINE(i,j) + f*bicubic(x,y,IMAGE)-sky
              
           case(3)
              ! bspline interpolation
              KOMBINE(i,j) = KOMBINE(i,j) + f*bspline_func(x,y,IMAGE)-sky
              
           end select

        enddo
     enddo
     if( sint == 3 ) then
        call bspline_shutdown
     endif

     if( onedate /= '' .and. oneexp > 0.0 ) then
        write(buf,'(a,1x,f0.2,1x,a,1x,f0.3,a,f0.3,1x,f0.2)',iostat=i) &
                   trim(onedate),oneexp,trim(filtr),f,'+-',fs,sky
        if( i /= 0 ) buf = trim(onedate)//' ..'
        etime = etime + oneexp
        if( first == '' ) first = onedate
        
        nalldate = nalldate + 1
        if( nalldate > maxalldate ) then
           allocate(alld(maxalldate))
           alld = alldate
           deallocate(alldate)
           maxalldate = 2*maxalldate
           allocate(alldate(maxalldate))
           alldate(1:maxalldate/2) = alld
           deallocate(alld)
        endif
        alldate(nalldate) = buf 
        
     endif
     n = n + 1
665  continue
     deallocate(IMAGE)

  end if

  enddo
20 continue

  if( n == 0 ) stop 'No input image.'
  if( bitpix == 0 ) bitpix = bpix
  bitpix = -32

  ! normalize
  select case (norm)
  case(1)
     KOMBINE = KOMBINE/n
  case(2)
     allocate(obr(nx*ny))
     obr = pack(KOMBINE,.true.)
     call rmean(obr,sky,y)
     KOMBINE = KOMBINE / sky
     deallocate(obr)
  end select

  ! range checking
  if( bitpix > 0 )then
     maxvalue = 2**bitpix - 1
     where( KOMBINE < minvalue )
        KOMBINE = minvalue
     endwhere
     where( KOMBINE > maxvalue )
        KOMBINE = maxvalue
     endwhere
  end if

! form the output image
  Inquire(file=namekombi,exist=exist)
  if( exist ) then
     open(26,file=namekombi,iostat=i)
     close(26,status='DELETE', iostat=i)
     if( i /= 0 ) then
        write(*,*)"Error occured during remove of output file."
        goto 666
     endif
  endif
  Call ftinit(26,namekombi,1,istat)
  naxis = 2
  naxes = (/nx, ny/)
  Call ftphps(26,bitpix,naxis,naxes,istat)
  if( bitpix > 0 ) then
     call ftpkye(26,'BSCALE',1.0,10,'',istat)
     call ftpkye(26,'BZERO',2.0**(bitpix-1),10,'',istat)
  endif
  Call ftpkys(26,'OBJECT',obj,'',istat)
  Call ftpkys(26,'FILTER',filtr,'',istat)
  call ftpkye(26,'EXPTIME',etime,10,'Sum of all times',istat)
  if( first /= '' ) call ftpkys(26,'DATE-OBS',first,'UTC of serie exposure start',istat)
  if( nalldate > 0 ) then     
     write(buf,'(a,i0,a)',iostat=i) &
          'The image is a result of composition of ',nalldate,' exposure(s).'
     if( i== 0 ) call ftpcom(26,buf,istat) 
     call ftpcom(26,'The time of start, the exposure time, the filter,',istat) 
     call ftpcom(26,'the intensity multiplier and relative sky value for each image used:',istat)
     do i = 1, nalldate
        call ftpcom(26,alldate(i),istat)
     enddo
     call ftpcom(26,'All times are derived from DATE-OBS and EXPTIME keywords.',istat)
  endif
  if( norm > 0 ) then
     select case (norm)
     case(1)
        buf = 'Fluxes are normalized by number of images.'
        i = 0
     case(2)
        write(buf,'(a,f0.3)',iostat=i)'Fluxes are normalized by sky value: ',sky
     end select
     if( i == 0 ) call ftpcom(26,buf,istat)
  end if
  select case (sint)
  case(0)
     buf = 'No interpolation has been used.'
  case(1)
     buf = 'Bilinear interpolation has been used.'
  case(2)
     buf = 'Bicubic interpolation has been used.'
  case(3)
     buf = 'Bspline interpolation has been used.'
  end select
  call ftpcom(26,buf,istat)
  call ftpcom(26,version,istat)
  if( istat /= 0 ) goto 666
  
  Call ftp2de(26,1,nx,nx,ny,KOMBINE,istat)
  
  CALL ftclos(26,istat)
  if( istat /= 0 ) goto 666

  ! print output info
  write(*,*) 'Output image: ',trim(namekombi)
  write(*,*) 'Dimension: ',nx,'x',ny

666 continue
  if( allocated(KOMBINE) ) deallocate(KOMBINE)

  if( istat /= 0 )then
     write(*,*) 'The following error(s) has been occured:'
     do
        Call ftgmsg(buf)
        if( buf == ' ' ) exit
        write(*,*) trim(buf)
     enddo
  endif

end program skombine

